#include <cblas.h> // 对于BLAS
#include <lapacke.h> // 对于LAPACK
#include <iostream>
#include <fstream>
#include <stdexcept>
#include <cmath> 
#include <vector>
#include "Sphere.hpp"
#include "Spline.hpp"
#include "U.hpp"

std::string directory = "src/build/";  // 目标目录

double pi = M_PI;

// 定义心形曲线的方程
double heartShape_x1(double t) {
    return sqrt(3) * cos(t);
}

double heartShape_y1(double t) {
    return 2.0 * sqrt(3) * sin(t) / 3 + 2.0 * pow(3, 0.25) / 3 * sqrt(std::fabs(cos(t)));
}

double heartShape_x2(double t) {
    return sin(t) + t*cos(t);
}

double heartShape_y2(double t) {
    return cos(t) - t*sin(t);
}

double heartShape_x3(double t) {
    return sin(cos(t))*cos(sin(t));
}

double heartShape_y3(double t) {
    return sin(cos(t))*sin(sin(t));
}

double heartShape_z3(double t) {
    return cos(cos(t));
}

std::vector<double> generateT(int m, int n){
    std::vector<double> t;
    if(n==1){
        double delta = 2.0 * M_PI / m; 
        for (int i = 0; i <= m; ++i) {
            t.push_back(i * delta + M_PI / 2);
        }
    }
    else if(n==2){
        double delta = 6.0 * M_PI / m; 
        for (int i = 0; i <= m; ++i) {
            t.push_back(i * delta);
        }
    }
    else if(n==3){
        double delta = 2.0 * M_PI / m; 
        for (int i = 0; i <= m; ++i) {
            t.push_back(i * delta);
        }
    }

    return t;
}

// 生成控制点
std::vector<std::vector<double>> generateControlPoints(int m, int n) {
    std::vector<std::vector<double>> controlPoints;
    if(n==1){
        double delta = 2.0 * M_PI / m; 
        for (int i = 0; i <= m; ++i) {
            double t = i * delta + M_PI / 2;
            double x = heartShape_x1(t);
            double y = heartShape_y1(t);
            controlPoints.push_back({x, y});
            /*if(m==10){
                std::cout << "x=" << x << ", y=" << y << std::endl;
            }*/
        }
    }
    else if(n==2){
        double delta = 6.0 * M_PI / m; 
        for (int i = 0; i <= m; ++i) {
            double t = i * delta;
            double x = heartShape_x2(t);
            double y = heartShape_y2(t);
            controlPoints.push_back({x, y});
            /*if(m==10){
                std::cout << "x=" << x << ", y=" << y << std::endl;
            }*/
        }
    }
    else if(n==3){
        double delta = 2.0 * M_PI / m; 
        for (int i = 0; i <= m; ++i) {
            double t = i * delta;
            double x = heartShape_x3(t);
            double y = heartShape_y3(t);
            double z = heartShape_z3(t);
            controlPoints.push_back({x, y, z});
        }
    }

    return controlPoints;
}

template <typename T>

double polynomial(double x, T p) {
    return p.SplineReturn(x);  
}

int main()
{
    // cumulative
    for (int i = 1; i < 4; i++){
        std::vector<int> array = {10, 40, 160};
        for (int m : array) {
            std::ofstream file(directory + "data_E" + std::to_string(i) + "_" + std::to_string(m) + "_cumulative" + ".csv");
            if (!file.is_open()) {
                std::cerr << "Failed to open file" << std::endl;
                return 1;
            }

            std::vector<std::vector<double>> controlPoints = generateControlPoints(m,i);

            std::vector<double> x1_t;
            std::vector<double> y1_t;
            std::vector<double> z1_t;
            std::vector<double> t1;

            if(i==1||i==2){
                std::vector<double> x1, y1;
                for (const auto& point : controlPoints) {
                    x1.push_back(point[0]);
                    y1.push_back(point[1]);
                }
                calculateCumulativeChordLengths(x1, y1, t1, x1_t, y1_t);

                if(i==1){
                    BSpline px(t1,x1_t,3,"periodic");
                    BSpline py(t1,y1_t,3,"periodic");

                    for (double t = t1.front(); t <= t1.back(); t += 0.002) {
                        double x = polynomial<BSpline>(t,px);
                        double y = polynomial<BSpline>(t,py);
                        file << x << "," << y << std::endl;
                    }
                }
                else{
                    BSpline px(t1,x1_t,4);
                    BSpline py(t1,y1_t,4);

                    for (double t = t1.front(); t <= t1.back(); t += 0.002) {
                        double x = polynomial<BSpline>(t,px);
                        double y = polynomial<BSpline>(t,py);
                        file << x << "," << y << std::endl;
                    }
                }
            }
            else if(i==3){
                std::vector<std::vector<double>> points_t;
                calculateCumulativeChordLengths(controlPoints, t1, points_t);
                for (std::vector<double> t : points_t) {
                    x1_t.push_back(t[0]);
                    y1_t.push_back(t[1]);
                    z1_t.push_back(t[2]);
                }

                BSpline px(t1,x1_t,3,"periodic");
                BSpline py(t1,y1_t,3,"periodic");
                BSpline pz(t1,z1_t,3,"periodic");

                for (double t = t1.front(); t <= t1.back(); t += 0.002) {
                    double x = polynomial<BSpline>(t,px);
                    double y = polynomial<BSpline>(t,py);
                    double z = polynomial<BSpline>(t,pz);
                    file << x << "," << y << "," << z << std::endl;
                }
            }

            file.close();
        }

        std::ofstream file(directory + "data_E" + std::to_string(i) + "_cumulative" + ".csv");
            if (!file.is_open()) {
                std::cerr << "Failed to open file" << std::endl;
                return 1;
            }
            std::vector<std::vector<double>> matrix;
            matrix = generateControlPoints(10000, i);
            if(i==1||i==2){
                for (std::vector<double> d : matrix) {
                    file << d[0] << "," << d[1] << std::endl;
                }
            }
            else if(i==3){
                for (std::vector<double> d : matrix) {
                    file << d[0] << "," << d[1] << "," << d[2] << std::endl;
                }
            }

            file.close();
    }

    // equidistant
    for (int i = 1; i < 4; i++){
        std::vector<int> array = {10, 40, 160};
        for (int m : array) {
            std::ofstream file(directory + "data_E" + std::to_string(i) + "_" + std::to_string(m) + "_equidistant" + ".csv");
            if (!file.is_open()) {
                std::cerr << "Failed to open file" << std::endl;
                return 1;
            }

            std::vector<std::vector<double>> controlPoints = generateControlPoints(m,i);

            std::vector<double> t1;

            if(i==1||i==2){
                std::vector<double> x1, y1;
                for (const auto& point : controlPoints) {
                    x1.push_back(point[0]);
                    y1.push_back(point[1]);
                }
                t1 = generateT(m,i);

                if(i==1){
                    BSpline px(t1,x1,3,"periodic");
                    BSpline py(t1,y1,3,"periodic");

                    for (double t = t1.front(); t <= t1.back(); t += 0.002) {
                        double x = polynomial<BSpline>(t,px);
                        double y = polynomial<BSpline>(t,py);
                        file << x << "," << y << std::endl;
                    }
                }
                else{
                    BSpline px(t1,x1,4);
                    BSpline py(t1,y1,4);
                    for (double t = t1.front(); t <= t1.back(); t += 0.002) {
                        double x = polynomial<BSpline>(t,px);
                        double y = polynomial<BSpline>(t,py);
                        file << x << "," << y << std::endl;
                    }
                }
            }
            else if(i==3){
                std::vector<double> x1, y1, z1;
                std::vector<std::vector<double>> points_t = controlPoints;
                t1 = generateT(m,i);
                for (std::vector<double> t : points_t) {
                    x1.push_back(t[0]);
                    y1.push_back(t[1]);
                    z1.push_back(t[2]);
                }

                BSpline px(t1,x1,3,"periodic");
                BSpline py(t1,y1,3,"periodic");
                BSpline pz(t1,z1,3,"periodic");

                for (double t = t1.front(); t <= t1.back(); t += 0.002) {
                    double x = polynomial<BSpline>(t,px);
                    double y = polynomial<BSpline>(t,py);
                    double z = polynomial<BSpline>(t,pz);
                    file << x << "," << y << "," << z << std::endl;
                }
            }

            file.close();
        }

        std::ofstream file(directory + "data_E" + std::to_string(i) + "_equidistant" + ".csv");
            if (!file.is_open()) {
                std::cerr << "Failed to open file" << std::endl;
                return 1;
            }
            std::vector<std::vector<double>> matrix;
            matrix = generateControlPoints(10000, i);
            if(i==1||i==2){
                for (std::vector<double> d : matrix) {
                    file << d[0] << "," << d[1] << std::endl;
                }
            }
            else if(i==3){
                for (std::vector<double> d : matrix) {
                    file << d[0] << "," << d[1] << "," << d[2] << std::endl;
                }
            }

            file.close();
    }

    return 0; 
}